function groups = findgroups( img, join, treshold );
% groups=findgroups( img, join=1, treshold = 0.1 ) - finds groups in image, if join, joins
% groups
% image is 2d matrix of vaslues in [0,1) .
% returns vector of structures with fields:
% centerX, centerY - group center of mass
% maxX, minX, maxY, minY - group boundaries
% weight - total value of pixels in group

% quck algorithm, suboptimal

% init
groups = [];
groupsnum = 0;
joins = []; % to column matrix - which groups should be joined

if exist('join') ~= 1
	join=1;
end

if exist( 'treshold' ) ~= 1
	treshold = 0.1;
end

% image size
w = size( img, 2 );
h = size( img, 1 );

% map of group, contains goup index for each pixel, or 0
map = zeros( h, w );

% scan image
for y=1:h
	for x=1:w
		value = img( y, x);
		% is current pixel active?
		if value >= treshold
			% add to group. find group index
			%newindex = size( groups, 1) + 1;
			indexleft = 0;
			indexabove = 0;
			% check left pixel
			if x > 1 && map( y, x-1) ~= 0
				indexleft = map( y, x-1);
			end
			% check pixels above
			if y > 1
				if map( y-1, x) ~= 0
					indexabove = map( y-1, x);
				elseif x > 1 && map( y-1, x-1) ~= 0
					indexabove = map( y-1, x-1);
				elseif x < w && map( y-1, x+1) ~= 0
					indexabove = map( y-1, x+1);
				end
			end
			
			% select index for the pixel
			index = 0;
			if indexabove > 0 
				index = indexabove;
				if indexleft > 0 && indexabove ~= indexleft
					%joing groups indexabove and indexleft
					joins = [ joins ; indexabove indexleft ];
				end
			elseif indexleft > 0
				index = indexleft;
			else
				% TODO create new group
				groupsnum = groupsnum + 1;
				index = groupsnum;
			end
				
			% apply index
			map( y, x ) = index;
			
		end % if above treshold
	end % x
end % y

% read groups from as struct with fields:
% - sumX, sumY - sum of X and Y coordinates
% -maxX, maxY, minX, minY
% - sumValue - sum of pixel value

nullsum = struct('sumX', 0, 'sumY', 0, 'sumValue',0, 'maxX', 0, 'maxY', 0, 'minY', 9999, 'minX', 9999 );
sums =repmat( nullsum, groupsnum, 1);

% scan map
for y=1:h
	for x=1:w
		groupid = map(y,x);
		value = img( y, x );
		if groupid ~= 0
			sum = sums(groupid);
			sum.sumX = sum.sumX + x*value;
			sum.sumY = sum.sumY + y*value;
			sum.sumValue = sum.sumValue + value;
			sum.maxX = max( x, sum.maxX );
			sum.maxY = max( y, sum.maxY );
			sum.minX = min( x, sum.minX );
			sum.minY = min( y, sum.minY );

			sums(groupid) = sum;
		end % if group
	end % x
end %


% join groups
if join ~= 0
	for j = 1:size(joins, 1 )
		g1 = joins(j,1);
		g2 = joins(j,2);
		
		% add sums
		sums(g1).sumX = sums(g1).sumX + sums(g2).sumX;
		sums(g1).sumY = sums(g1).sumY + sums(g2).sumY;
		sums(g1).sumValue = sums(g1).sumValue + sums(g2).sumValue;
		
		% get extremes
		sums(g1).maxX = max( sums(g2).maxX, sums(g1).maxX );
		sums(g1).maxY = max( sums(g2).maxY, sums(g1).maxY );
		sums(g1).minX = min( sums(g2).minX, sums(g1).minX );
		sums(g1).minY = min( sums(g2).minY, sums(g1).minY );
	
		% reset the second group
		sums(g2) = nullsum;
	end
end %join

% read groups from map
for i = 1:size(sums,1)
	sum = sums(i);
	% convert to group if not null
	if sum.maxX ~= nullsum.maxX
		group = sum; % convienient - group extends sum
		group.centerX = sum.sumX / sum.sumValue;
		group.centerY = sum.sumY / sum.sumValue;
		
		groups = [ groups ; group ];
	end
end % i


					
				
				
			
				
			
